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Cross  Validated  Spline  Methods  for  Direct  and  Indirect  Sensing  Experiments 

Grace  Wahba 

University  of  Wlsconsin-Madison 
Statistics  Department 

Abstract 


We  first  give  an  overview  of  cross-validated  smoothing  spline  methods 
for  the  analysis  of  data  from  direct  and  indirect  sensing  experiments. 

Then  we  summarize  some  results  involving  two  particular  applications  of 
the  approach  -  (1)  a  deconvolution  problem  with  nonnegativity 
constraints  on  the  solution;  (2)  estimation  of  horizontal  divergence  and 
vorticity  from  discrete,  scattered,  global  measurements  of  the  upper  air  wind 
field. 


1. 


For  several  years,  we  have  been  developing  data  analysis  techniques  for  a 
fairly  general  class  of  (direct  and  Indirect)  sensing  experiments.  The  model 
may  be  described  as  follows:  f  Is  seme  unknown  function  of  one  or  more 
variables,  which  Is  assumed  to  be  in  a  (Hilb«rt)  space  H  of  "smooth"  functions, 
and  one  observes  {y^} 

y^  *  Ljf  +  e<  1  «  l,2,...fn,  (1.1) 


Overview  of  Cross  Validated  Spline  Methods  for  Direct  and  Indirect 
Sensing  Experiments. 


where  the  data  functionals  L.,,...,L_  are  n  bounded  linear  fi'nctlonals  on  H 

i  n 

and  the  are  Independent,  zero  mean  measurement  errors,  with  variances 
Wja2,  1  *  l,2,...,n.  The  parameter  a2  may  be  unknown.  Examples  of  are 
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Example  (1)  corresponds  to  the  "direct  sensing"  problem,  example  (2)  to  a 
Fredholm  Integral  equation  of  the  first  kind  given  discrete,  noisy  data,  and 
example  (3)  corresponds  to  a  partial  differential  equation.  Example  (2)  Is 
pervasive  In  many  problems  In  geophysics,  meteorology  and  medicine. 

Usually  one  desires  to  estimate  f(P),  for  all  Pe  some  set  S.  Sometimes  only 
moments  /w^(P)f(P)dP  are  desired.  In  other  applications  derivatives  of  f  may  be 
desired.  Side  Information  may  freq!?air.tly  be  available,  such  as 

0  <  f(P)  positivity 

a(P)  <  f(P)  <  b(P)  boundedness  ^  2) 

0  <f'(P)  monotonicity 

0  <f"(P)  convexity 

etc.  More  generally  we  have  considered  the  case  where  feCCH,  where  C  Is  any 
closed  convex  set  In  H  which  Is  the  Intersection  of  a  "smoothly  varying"  family  of 
hyperplanes.  See  [2],  H  can  be  chosen  so  that  there  exists  a  C  with  the 
requisite  properties  for  each  of  the  examples  In  (1.2). 

The  general  approach  to  the  estimation  of  f  we  have  developed  and  tested 
goes  as  follows:  The  estimate,  call  It  fy  Is  the  solution  to  the  minimization 
problem:  Find  feH  to  minimize 

<’-3> 

subject  to  feC.  Here  J  .  Is  a  sewlnorm  on  H,  or  "penalty  function" 

m  *o 

Indexed  by  the  parameters  m  and  6,  (to  be  described  In  a  monent),  and  X  Is  the 
bandwidth  or  smoothing  parameter.  Some  examples  of  J  are 


CD  JL(f)  -  /(flB,(P))*dP 


(2)  J2  3(f)  •  77(fxx*+26fxy*+6*fyy*)dxdy  or,  In  d  dimensions 


I  j-  a  i  /-•/(—*-- 


),dxr.-.dxd. 


(3)  J.Cf) 


J  (^f)2dP  ' 

sphere 


where  A  Is  the  Laplace-Beltraml  operator  on  the  sphere 

Af  *  — - —  f..  ♦  —rrCcos^f.).  ,  X  *  longitude,  +  *  latitude. 

Co$*^  AA  co*t  T  f  ' 

Another  example  appears  In  Section  3  below,  where  the  parameter  6  plays  a 
different  role.  The  Important  bandwidth  parameter  X,  and  sometimes  m  and  6', 
can  be  estimated  by  the  method  of  generalized  cross  validation  (GCV).  The 
6CV  method  has  recently . been  extended  to  cover  problems  where  the' 
constraints  feC  are  Imposed  [21,24].  In  example  (1),  If  the  L^f  are  point 
evaluation  functtonals,  that  Is',  L^f  *  f(P.|),  then  f^  Is  a  polynomial  smoothing 
spline*  In  example’ (2)  f^  will  be  a  thin  plate  spline  [11,13]  and  In  example 
(3)  f,  will  be  a  spline  on  the  sphere  [20]. 

f^  can  be  shown  to  be  a  Bayes  estimate,  with  a  certain  prior  determined 
by  H  and  J  ([1,6]).  On  the  other  hand.  It  can  be  shown  that  f^  Is  a  natural 
generalization  of  the  output  of  a  low  pass  filter,  see  [7].  If  the  t^f’s 
are  derivatives,  and  the  limit  Is  taken  as  X*0,  then  spline  collocation  methods 
for  the  solution  to  differential  equations  results  [10]. 

Table  1  gives  a  list  of  some  of  the  problems  Involved  In  estimating  f  by 
cross  validated  spline  methods,  and  relevant  references  of  the  author,  her 

A 

collaborators,  and  students.  In  Table  1,  X  Is  the  GCV  estimate  of  X. 


ft 

iii 


TabU  1 


\  •  \ 

Problems  In  the  estimation  of  f ,  with  references 

•  \  ' 

(1)  Explicit  representation  of  the  minimi zer  of  (1.3)  In  various  contexts 
0  *4,5,11 ,13,18,20,22]. 

(2)  Generalized  cross  validation  methods  for  choosing  Xjn*  and  6  [4,7,8,21 ,24]. 

(3)  Use  of  prior  Information  to  choose  H  and  J  [22,23]. 

(4)  Smoothing  of  multidimensional  scattered  data  [11,13]. 

(5)  Smoothing  splines  on  the  sphere  [20,22]. 

(6)  Numerical  methods-for  computing  f£  with  large  data  sets  [14,21,25]. 

(7)  Quadrature  methods  tailored  to  specific  Integral  equation  problems  [19,28]. 

(8)  Optimal  design  (choice  of  tj,...,tn)£l0,17],  - 

(9)  Convergence  properties  of  fj  [3,4,7,12]. 

(10)  Convergence  properties  when  the  constraints  are  discretized  [2]. 

(11)  Confidence  Intervals  for  fj  [26]. 

(12)  Relation  of  fj  to  Bayes  estimates  and  Weiner  filtering  [1,6]. 

(13)  Nonlinear  data  functionals  [29]. 

*  > 

Four  expeclally  Interesting  areas  of  application  are 
.  (1)  Smoothing  of  scattered,  noisy  data  In  two  and  higher  dimensions 
with  thin  p.1ate  splines  [11,13].  Transportable  code  Is  available 
from  the  Madison,  WI  Academic  Computing  Center. 

(2)  Approximate  solution  of  Abel's  Integral  Equations,  as  they  occur  In 

stereology  and  computerized  tomography  [16,28].  These  equations  have 
a  singular  kernel.  r' 

(3)  Equations  of  radiative  transfer  as  they  occur,  for  example  In  the 
estimation  of  vertical  temperature  profiles  from  satellite-observed 
radiances  [29].  These  equations  are  mildly  nonlinear. 

(4)  As  an  alternative  to  logistic  regression  [27]. 


Tilt  Integral  equation  Methods  have  found  use  In  a  number  of  engineering 
applications,  for  example  In  the  determination  of  adsorption  energy 
distributions  from  an  adsorption  Isotherm  experiment  [30]  and  In  the  recovery 
of  aerosol  size  distributions  from  Harple  Impactor  data  [31]. 

For  lack  of  space  we  have  omitted  explicit  mention  of  recent  related 

f" 

work  by  others •  In  particular  by  Cox;  Chow,  Geman  and  Wu;  Rosenblatt; 

L11;  Silverman;  Speckmanj and  Utreras. 

In  the  remainder  of  this  paper  we  briefly  review  recent  nunerlcal  results 
In  two  Interesting  areas:  Section  2:  Deconvolution  with  positivity  constraints 
and  Section  3:  Estimation  of  wind  horizontal  divergence  and  vortlclty  from 
scattered  noisy  wind  vector  measurements. 


2.  Deconvolution  with  Nonneoatlvlty  Constraints 
He  numerically  studied  the  convolution  equation 


1  , 

y1  *  ^(j;-$)f(s)ds  ♦  ei,  .  1  *  l,2,...,i 


(2.1) 


f(s)  >  0,.0  <  s  <  1 


I 

with  J(f)  «  ^(f"(s))2ds.  The 'results  appear  In  [24].  The  cross  validated  spline 
method  prescribes  the  estimate  of  the  solution,  call  It  f^,  say,  as  the  ninlmlzer 

Mi  +  xl<f"(s»,ds  (*•*> 


In  the  nonnegative  quadrant  of  a  certain  Sobolev  space  M22,  (which  we  will 
not  discuss  further). 

For  n  of  moderate  size  (say  32-256)  a  good,  readily  computable  approximation 
to  the  mlnlmlzer  f ^  of  (2.2)  may  be  found  by  approximating  f^  by  a 
trlgnometrlc  polynomial  of  the  form 

n/2  n/2-1 

f.(x)  *  cu  ♦  I  olCmZtm  ♦  l  6 1  slnZmm 
a  v»1  v-1 

and  finding  the  a's  and  0's  to  minimize  (2.2)  subject  to  the  discretized 
positivity  constraints 
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f(jj!  >  Oj  1*  ■  1,2,...,n. 

For  fixed  X,  a  quadratic  programming  problem  with  n  unknowns  and -n  linear.  - 

Inequality  constraints  results.  GCV  for  constrained  problems  may  be used 

•ip 

to  choose  A;  for  each  trfalvalueofA*  q.p.  must  be  solved,  a  good  starting 
guess  can  be  obtained  by  applying  GCV  to  the  unconstrained  problem.  (Thus, 
thermethodis  not  "cheaper  however  It' can  be  very  cheap  compared  tO'-the 
cost  of  some  experiments.!. 

He  give  here  a  single  numerical  example,  from  [24].  Fig.  2.1  gives  a 

* 

plot  of  the  convolution  kernel  k(t)  (assuned  periodic).  Fig.  2.2  gives 
a  plot  of  the  (true)  test  solution  f(x),  0  <  x  <  1,  "exact"  data  g(x), 

1 

g(x)  ■  ^k(x-s)f(s)ds  0  <  x  <  1 

and  simulated  measured  data 

> 

yi *  *(»>  +  *r  1 " 1,2 . n 

for  n  ■  64.  The  were  Independent  normally  distributed  pseudorandom 
varlables-wlth  cannon  standard  deviation  o  *  .05.  The  fact  that  there  are 
two  distinct  peaks  In  the  true  f  Is  not  at  all  obvious  from  the  measured  data. 
Fig.  2.3  gives  the  true  f  (again)  and  It  gives  f^,  which  is  the  crossvall dated 
spline  estimate  of  f  obtained  without  Imposing  positivity  constraints.  Fig.  2.3 

f  * 

also  gives  f?  ,  which  Is  the  estimate  for  f  with  the  nonnegativity  constraints 
imposed.  Ag  Is  the  GCV  estimate  of  A  for  constrained  problems. 

The  unconstrained  cross-va 11  dated  spline  estimate  Is  not  bad,  however, 
spurious  oscillations  are  clearly  visible.  The  constrained  solution  not 
only  eliminates  the  spurious  oscillations.  It  enhances  the  resolution  of  the 
two  peaks.  This  Innocuous  looking  problem  Is  lllposed  to  a  high  degree. 

V 
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If  there  were  no  measurement  error,  there  would  be  only  about  20  linearly 
independent  pieces  of  information  in  the  64  dimensional  data  vector  recorded 
to  8  significant  figures.  The  nonnegativity  constraints  are  adding  important 
information.  For  further  Information  concerning  the  relative  degree  of 
lllposedness,  see  [21].  Unsophisticated  estimation  methods  will  give  garbage. 

3.  Vector  Fields  on  The  Sphere 

In  meteorology  and  geophysics  measurements  of  vector  fields  are  made 
at  discrete  locations  around  the  world  and  it  is  desired  to  estimate  the 
global  vector  field.  In  meteorology,  for  example,  the  wind  field  at  the 
500  millibar  height  (height  at  which  the  pressure  Is  500  millibars)  and  other 
heights,  is  measured  at  the  worldwide  radiosonde  network.  It  is  desired  to 
estimate  the  vortlclty  (curl)  c  and  the  divergence  D  of  this  wind  field  on 
a  regular  grid,  to  be  used  as  initial  conditions  for  the  differential 
equations  of  numerical  weather  forecasting.  The  following  is  abstracted 
from  [22]. 

Letting  U(P)  and  V(P)  be  the  east  and  north  wind  respectively  at 
P  *  (A,<J>),  the  data  are 

Ui  -  U(P1)  +  c”,  Vi  *  V(P1)  +  e J,  (3.1) 

U  V 

where  and  ej  are  measurement  errors,  and  it  is  desired  to  estimate  c  and 
0  on  a  regular  grid,  ;  and  0  are  related  to  U  and  V  by 

«  *  SHjt-  ♦  ft  "  ’  f  +  <3'2> 


where  a  is  the  radius  of  the  earth.  There  exists  (by  Helmoltz  Theorem) 
two  functions  'F(P)  and  *(P),  called  the  stream  function  and  the  velocity 
potential  respectively,  with  the  following  properties: 


r.eaw'a. 


Me  estimate  Z  andOf roar  Ilf  and  V1  by  solving  a  minimization 
problem  of  the  fore:  Find  f  and  P  In  an  appropriate  function  space  to 
'minimize  ~  .  - — 

♦  X[J0,(t>  ♦  jo(2)(«)] 


(3.5) 


and  then  obtaining  c  and  D  analytically  as  A?  and  Ad.  The  parameter  \  Is 
the  usual  bandwidth  parameter  and  the  parameter  6  controls  the  relative 
amount  of  energy  In  the  divergent  and  nondlvergent  part  of  the  wind.  An 
approximation  to  the  mlnlmlzer  of  (3.5)  may  be  obtained  by  first  approximating 
¥  and  «as  a  finite  number  of  spherical  harmonics 

M  A  .  HI 

f  *  A  2  & sYt*  ♦  ■  1,  £  ®*sV 

*  ,  AfI  $»-£  **  *  £-1  S*-£ 


(3.6) 


The  spherical  harmonics  (Y*)  are  the  eigenfunctions  of  the  Laplace -Beltrami 
operator  A, 


1  AyJ  ■  -  A(i+l)Yj, 

and  approximating  ¥  and  P  this  way  Is  analogous  to  the  approximation  of 

p 

f  In  Section  2  by  a  finite  number  of  sines  and  cosines.  It  can  be  shown  that 
|  If  and  are  Isotropic  seminorms  (Isotropic  ■  unchanged  under  arbitrary 

I  rotations  of  the  coordinate  system)  then  for  ¥  and  P  of  (3.6) , 


(3.7) 


for  som  {Xjjy,X^}.  For  sore  details,  and  in  partlcular.for  a  discussion 
concerning  the  choice  of  thef^}  fron  historical  data,  see  [22]. 

By  substituting  (3.7)  Into  (3.5),  for  fixed  A,6,  the  Minimization 
problem  reduces  to  finding  the  {a^ ,Bts )  which  minimize  a  quadratic  font. 

The  parameters  A  and  6  are  estimated  by  6CV. 

/  To  see  how  wefL.thls  method  may  be  expected  to  do  on  real  neteorolcglcal 
data,  realistic  "true"  500  millibar  stream  function-velocity  potential 
pairs  (f,4)  were. generated  and  the  "true"  wind  fields  at  114  North  American 
weather  stations  determined  using  (3.3).  Realistic  measurement  errors.. 

(s.d.  In  each  component  of  2.5  meters/sec)  were  added.  Fig.  3.1  gives  the 
simulated  wind  data,  and  Fig.  3.2  gives  the  estimated  wind  field,  which  is 
obtained  by  using  (3.3)  In  conjunction  with  the  estimated  (Y.ft).  .Figs. 

3.3  and  3.4  give  the  "true"  and  estimated  vortlcity  and  divergence 

respectively.  The  results  appear  to  be  excellent  when  compared  to  previous 

I 

estimation  methods.  One  of  the  results  of  this  experiment  was  that  the 

•  !  - 

estimates  were  sensitive  to  changes  In  S  as  well  as  A  and,  the  GCV  methcd 
gives  a  good  estimate  of  an  optimal  6  as  well  as  A.  In  a  problem  like  this. 
It  can  be  very  Important  to  parameterize  J  In  a  physically  meaningful  way, 
as  well  as  In  a  mathematically  well  posed  May.  Methods  for  doing  this  are 
discussed  In  [22],  see  also  [15]. 
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